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Abstract 

Background: Gen(om)e duplication events are hypothesized as key mechanisms underlying the origin of 
phenotypic diversity and evolutionary innovation. The diverse and species-rich lineage of teleost fishes is a renowned 
example of this scenario, because of the fish-specific genome duplication. Gene families, generated by this and other 
gene duplication events, have been previously found to play a role in the evolution and development of innovations 
in cichlid fishes - a prime model system to study the genetic basis of rapid speciation, adaptation and evolutionary 
innovation. The distal-less homeobox genes are particularly interesting candidate genes for evolutionary novelties, such 
as the pharyngeal jaw apparatus and the anal fin egg-spots. Here we study the dlx repertoire in 23 East African cichlid 
fishes to determine the rate of evolution and the signatures of selection pressure. 

Results: Four intact dlx clusters were retrieved from cichlid draft genomes. Phylogenetic analyses of these eight dlx loci 
in ten teleost species, followed by an in-depth analysis of 23 East African cichlid species, show that there is disparity in 
the rates of evolution of the dlx paralogs. Dlx3a and dlx4b are the fastest evolving dlx genes, while dlxla and dlx6a 
evolved more slowly. Subsequent analyses of the nonsynonymous-synonymous substitution rate ratios indicate that 
dlx3b, dlx4a and dlx5a evolved under purifying selection, while signs of positive selection were found for dlxla, dlxla, 
dlx3a and dlx4b. 

Conclusions: Our results indicate that the dlx repertoire of teleost fishes and cichlid fishes in particular, is shaped by 
differential selection pressures and rates of evolution after gene duplication. Although the divergence of the dlx 
paralogs are putative signs of new or altered functions, comparisons with available expression patterns indicate that 
the three dlx loci under strong purifying selection, dlx3b, dlx4a and dlx5a, are transcribed at high levels in the cichlids' 
pharyngeal jaw and anal fin. The dlx paralogs emerge as excellent candidate genes for the development of evolutionary 
innovations in cichlids, although further functional analyses are necessary to elucidate their respective contribution. 

Keywords: Distal-less homeobox gene, Molecular evolution, Cichlid fishes, Teleost fishes, Positive selection, Differential 
selection, Gene duplication, dN/dS 



Background 

Teleost fishes (Teleostei) are among the most diverse 
lineages on Earth and with nearly 30,000 species the most 
species-rich vertebrate group. This is in stark contrast 
to the more basal non-teleost ray-finned fishes that are 
characterized by small numbers of species. A causal explan- 
ation for this discrepancy in speciation rates between 
the derived Teleostei and the non-teleost ray-finned fishes 
might be the fish-specific genome duplication (FSGD) that 
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occurred in the ancestor of modern teleosts ([1-4] and ref- 
erences therein). It has been hypothesized that the FSGD 
has laid down the genetic conditions necessary for the 
evolution of phenotypic diversity [5], although the exact 
causes of diversification of such a large clade are likely to 
be more complex and most probably also include other 
factors [6]. 

The Hox gene clusters, which evolved through both tan- 
dem and whole genome duplications, represent illustra- 
tive examples for the contribution of duplicated genes 
to morphological evolution across the animal kingdom 
(see e.g., [7-9]). Together with other homeotic genes, 
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Hox genes play a crucial role in the development of 
the multicellular body plan (e.g., anterior-posterior pat- 
terning; [10]). Furthermore, Hox genes are known to be 
involved in the development of evolutionary novelties, 
such as walking limbs and the wings of insects [11-15]. It 
has been shown that different mechanisms such as cis- 
regulatory evolution, changes in protein function and post- 
transcriptional regulation of the Hox genes contribute to 
morphological diversification (reviewed in e.g., [8,15,16]). 

East African cichlid fishes show a remarkable level of 
phenotypic diversity between closely related species and 
constitute the most diverse adaptive radiations known 
[17-21]. Although several smaller radiations of cichlid fishes 
exist outside of Africa (e.g., in Central and South America), 
an astonishingly high number of cichlid species (close to 
1900 species [22]) evolved in and around lakes Malawi, 
Victoria and Tanganyika in the last few million to several 
thousand years [23,24]. The various cichlid species differ 
in body shape, coloration, reproductive biology and mouth 
morphology [25-27] - traits which are thought to, at 
least partly, underlie the evolutionary success of cichlid 
fishes [18,27,28]. Furthermore, several morphological inno- 
vations are unique to cichlids or specific lineages thereof. 
The highly modified and morphological diverse pharyngeal 
jaw apparatus, for example, correlates with the diversity in 
foraging strategies exploited by the different cichlid species 
[27-29]. The occurrence of several color morphs within 
species, sexual color dimorphism and anal fin egg-spots 
are three characteristic features of the extremely species- 
rich and mouthbrooding haplochromine lineage [30] . 

As a result of their great phenotypic diversity and high 
number of species, cichlid fishes provide an ideal set up 
to examine the genetic basis of rapid speciation, evolu- 
tionary innovations and adaptation [21,31-37]. An import- 
ant strategy is the study of so-called candidate genes, 
i.e., genes with known functions in development in other 
organisms such as zebrafish. For example, it has been 
shown that csflra, which was identified as xanthophore 
marker in zebrafish [38,39] is involved in the morpho- 
genesis of the egg-spots of haplochromine cichlids [31]. 
Furthermore, species-specific jaw shapes of different cich- 
lid species correlate with differences in early bmp4 ex- 
pression patterns, a gene which has also the potential to 
change the mandibular morphology in zebrafish [40]. 
Many of these candidate genes belong to larger gene 
families such as the endothelin family of ligands and re- 
ceptors that are putatively involved in the morphogen- 
esis of the pharyngeal jaw apparatus and pigmentation 
[36], and the above mentioned Hox gene clusters [41]. 

Recendy, Renz et al. [35] characterized seven distal-less 
homeobox (dlx) genes and examined their expression pat- 
terns in the developing pharyngeal arches and/or pharyngeal 
teeth of the haploclvomine cichlid Astatotilapia burtoni. 
The vertebrate dlx genes are widely known for their 



crucial roles in the development (of components) of the 
nervous system, craniofacial skeleton and connective tis- 
sue and in the formation of appendages [reviewed in 42]. 
These functions seem to be conserved across a wide range 
of animal taxa. For example, the vertebrate dlx genes are 
homologs of, and share several functions with, the single 
Distal-less (dll) gene of Drosophila [42]. Within verte- 
brates, the expression patterns of dlx homologs is similar 
in early development [35,42-45]. At the same time, dlx 
genes have been implicated with evolutionary novelties 
such as the eyespots in various butterfly species [46-48], 
the insect antenna [49,50] and the vertebrate craniofacial 
bones [51]. 

Phylogenetic analyses and the chromosomal arrange- 
ments of the vertebrate dlx genes suggest that the extant 
dlx repertoire has evolved by an initial tandem duplica- 
tion, followed by two rounds of whole genome duplica- 
tion in the lineage towards vertebrates and a third one in 
the lineage towards teleost fishes, the FSGD [35,44,45,52]. 
These duplication events resulted in multiple so-called dlx 
clusters, in which two dlx genes are located in a tail-to-tail 
arrangement on the respective chromosome. Linked dlx 
genes are transcribed coincidently due to shared de- 
regulating elements in the intergenic regions [35,42,43]. 
Four of these dlx clusters have been identified in teleost 
fish; dlxla-dlx2a, dlx3a-dlx4a, dlx3b-dlx4b and dlxSa- 
dlx6a [43,44]. Seven of these dlx genes have been identi- 
fied in the cichlid A. burtoni, where they are expressed in 
tissues that make up putative evolutionary innovations [35]. 

Here, we analyzed the dlx repertoire and diversity in 
detail in a phylogenetically representative set of 23 East 
African cichlid species in order to study the molecular 
evolution of this prominent developmental gene family. 
To this end, we first performed phylogenetic compari- 
sons of the dlx proteins including the sixty amino acids 
long homeobox domain in a range of teleost fishes in 
combination with blast searches of these sequences against 
the draft genomes of four cichlid species. Teleost and 
cichlid-specific phylogenies were examined to compare 
the rate of evolution of both between and within dlx 
gene trees. Several studies have shown that loci putatively 
involved in evolutionary innovations are characterized by 
adaptive protein evolution in cichlids [31,36,53]. There- 
fore, all loci were screened for elevated rates of protein 
evolution by means of d N /d s analyses. Our analyses in- 
dicate the presence of dlx3a in cichlids and that the dlx 
repertoire of cichlid fishes is shaped by differential selection 
pressures and rates of evolution, with signs of positive se- 
lection on specific sites in dlxla, dlx2a, dlx3a and dlx4b. 

Methods 

Dlx protein sequence comparison in teleost fishes 

The sequences of nine dlx proteins (i.e., dlxla, dlx2a, dlx2b, 
dlx3a, dlx3b, dlx4a, dlx4b, dlx5a and dlx6a) of seven teleost 



Diepeveen et al. BMC Evolutionary Biology 2013, 13:153 
http://www.biomedcentral.com/1471 -21 48/1 3/1 53 



Page 3 of 12 



species (i.e., zebrafish (Danio rerio), Atlantic cod (Gradus 
morhua), three-spined stickleback (Gasterosteus aculeatus), 
spotted green pufferfish (Tetraodon nigroviridis), Japanese 
pufferfish (Takifugu rubripes), Japanese medaka (Oryzias 
latipes) and Nile Tilapia (Oreochromis niloticus)) were 
obtained from Ensemble (release 68, July 2012; see 
Additional file 1 for accession numbers). Dlx2b was 
excluded from all further analyses, due to its lineage- 
specific loss in percomorphs, to which all studied species 
belong except D. rerio and G. morhua (see [35]). Se- 
quences were aligned with Tcoffee [54,55], ambiguous 
sites were removed and tblastx searches were performed 
to determine dlx protein sequences in the draft cichlid 
genomes of Astatotilapia burtoni, Neolamprologus brichardi 
and Pundamilia nyererei (BROAD Institute, unpublished 
data; see Additional file 1 for scaffold numbers). To deter- 
mine the rate of evolution for each of the dlx proteins, 
phylogenetic analyses were performed in PAUP* 4.0 [56] 
under parsimony settings and the number of amino acid 
changes was obtained. D. rerio or G. morhua was used as 
outgroup species and bootstrap analyses with 100 replicates 
were conducted to test the robustness of the obtained 
topologies. Next, the sixty amino acids long homeobox 
domain was extracted from the sequences and aligned 
to the homeobox domain of the single Distal-less (Dll) 
gene of Drosophila melanogaster [Ensemble: FBgn0000157] 
in Geneious 5.6 [57] for closer inspection of the con- 
servation of the domain and to identify gene-specific 
substitutions. 

Cichlid samples and genomic DNA sequencing 

White muscle and/or fin clip samples were collected 
during fieldwork in Zambia in 2007 and 2008 using a 
standard operating procedure described in [29]. In total 
23 Lake Tanganyikan cichlid species were included in this 
study (Additional file 2). Genomic DNA was extracted 
following a standard Proteinase K protocol [58]. Cichlid- 
specific PCR primers were designed based on available 
and/or draft genomic and transcriptomic cichlid sequences, 
which were identified by tblastx searches of publicly 
available dlx sequences from other teleost species (see 
Additional file 1 for species and accession numbers). This 
was done for eight dlx loci: dlxla, dlx2a, dlx3a, dlx3b, 
dlx4a, dlx4b, dlxSa and dlx6a (see Additional file 3 
for primer sequences). Standard PCR reactions, purifi- 
cation steps and sequencing reactions were set up and 
performed as described elsewhere [36]. PCR products of 
the partially sequenced loci were visualized with GelRed 
(Biotium) on a 1.5% agarose gel and sequenced on a 
3130x/ capillary sequencer (Applied Biosystems). Partial 
sequences were aligned and visually inspected using Codon 
Code Aligner 3.7.1 (CodonCode Corporation, Dedham, 
MA). Exon/intron boundaries were determined by hom- 
ology comparisons with the sequences from the other 



teleost species. All generated cichlid dlx sequences have 
been deposited into GenBank [GenBank: KC285366- 
KC285546] (Additional file 2). 

Phylogenetic analyses of cichlid samples 

Individual gene trees were constructed using maximum 
likelihood in PAUP* 4.0 [56] and Bayesian Inference in 
MrBayes 3.2 [59,60]. The best-fitting model of nucleotide 
substitution was determined with the corrected Akaike 
information criteria and likelihood ratio tests conducted 
in jModeltest 0.1.1 [61,62]. Bootstrap analyses with 100 
replicates were performed in PAUP* and MrBayes was 
run for 10.500.000 generations. Oreochromis tanganicae 
was used as outgroup (see e.g., [63]). Phylogenetic ana- 
lysis of a concatenated dataset of 9.2 kb was performed 
as described above in PAUP* to generate a common in- 
put tree file (treeBASE submission 14433) for the subse- 
quent analyses. 

Coding sequence data of the 23 cichlid species (treeBASE 
submission 14433) was assessed with both site- and branch- 
site models as implemented in the program Codeml of 
the software package PAML (Phylogenetic Analysis by 
Maximum Likelihood) 4.3 [64,65]. The following parame- 
ters were estimated for all eight dlx datasets under differ- 
ent models: the nonsynonymous/synonymous substitution 
rate ratio, co, the proportion of sites assigned to an w cat- 
egory, po,i,2> an( l the p and q parameters of the /3 distri- 
bution. Tests of positively selected sites were conducted 
by performing Likelihood Ratio Tests (LRT) of the fol- 
lowing model comparisons: Mia (Nearly Neutral) with 
M2a (Positive Selection), M7 (|3) with M8 (|3 & co s > 1), 
and M8a (p & co s = 1) with M8. The comparison between 
M0 (one-ratio) and M3 (discrete) was used as a test of 
variable b> among sites. The naive empirical Bayes (NEB; 
[66,67]) and the Bayes empirical Bayes (BEB; [68]) criteria 
were used to calculate the posterior probabilities for site 
classes and the BEB was used to identify sites under positive 
selection when the LRT was significant. To test whether 
the dlx genes evolved under non-neutral evolution in spe- 
cific lineages a LRT between the null model (co s = 1) and 
the alternative model (co s > 1) was performed in the branch- 
site analyses. Branches of interest, or so-called foreground 
branches, were chosen based on the results of the phylo- 
genetic analyses and branch tests performed in Hyphy 
([69], following [36]). 

Additional tests of positive selection on the partial dlx 
sequences were performed with the Sitewise Likelihood 
Ratio estimation of selection program (SLR; [70]) vl.3. 
The common input tree file was used (see above) and 
the significance level was set to 95%. 

Amino acid substitutions were screened for possible ef- 
fect on protein function with the program SIFT (Sorting 
Intolerant from Tolerant; [71]. 
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Results 

Dlx protein sequence comparison in teleost fishes 

The tblastx searches of the teleost dlx proteins resulted 
in the retrieval of eight dlx genes in all four cichlid spe- 
cies. Furthermore, the genomic locations of these dlx 
loci (Additional file 1) indicate that four dlx clusters are 
present in the cichlid lineage: dlxla-dlx2a; dlx3a-dlx4a; 
dlx3b-dlx4b and dlx5a-dlx6a. All other teleost species 
examined contain this full set of genes, except zebrafish, 
in which dlx3a could not be located, and medaka, in 
which dlx4b is missing, as previously noted [35,44,45]. 
Interestingly, in contrast to Renz et al. [35] we do find 
evidence for the existence of dlx3a in cichlids, including 
A. burtoni (Figure 1, Additional file 4). 

The sixty amino acid long homeobox domain of the 
eight teleost dlx proteins are highly conserved among 
teleost fish and even between teleosts and the single Dll 
protein of D. melanogaster (Additional file 4). Despite 
the high level of conservation, several locus-specific amino 
acid substitutions are present in the paralogs, making it 



possible to distinguish between individual dlx homeo- 
box domains. 

Phylogenetic analyses of the dlx protein sequences were 
performed to examine the rate of evolution of the dlx 
paralogs in teleost fishes. The overall and relative longest 
trees were found for dlx4b and dlx3a, while for dlxla and 
dlx6a the shortest tree lengths were observed (Figure 1 
and Table 1). Typically the longest branches were observed 
in the two basal species D. rerio and G. morhua. Interest- 
ingly, relatively long branch lengths for the branch towards 
the four cichlid species were observed for dlx3a and dlx6a, 
indicting elevated rates of molecular evolution. The oppos- 
ite scenario was observed in the overall more conserved 
dlxla and dlx5a proteins. To study these effects in more 
detail cichlid specific gene trees were constructed. 

The rate of dlx gene evolution in East African cichlid fishes 

To reconstruct the molecular evolutionary history of the 
dlx homologs in East African cichlid species, we deter- 
mined the rate of evolution and the signatures of selection 
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Figure 1 Maximum Likelihood phylogenetic hypotheses for the eight dlx paralogs in teleost fishes. (A) Dlxla (254 amino acids (aa)). 
(B) Dlx2a (276 aa). (C) Dlx3a (307 aa). (D) Dlx4a (259 aa). (E) Dlx3b (283 aa). (F) Dlx4b (257 aa). (G) Dlx5a (285 aa). (H) Dlx6a (247 aa). Bootstrap probabilities 
(PAUP*) above 50% are shown. 
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Table 1 Overall and relative tree lengths of teleost 
protein phylogenies 



dlx paralog 


Length 
(# amino acids) 


Tree length 
(absolute) 


Tree length 
(relative) 


dlxla 


254 


93 


0.366 


dlx2a 


276 


122 


0.442 


dlx3a 


307 


179 


0.583 


dlx3b 


283 


108 


0.382 


dlx4a 


259 


123 


0.475 


dlx4b 


257 


222 


0.864 


dlx5a 


285 


125 


0.439 


dlx6a 


2-1/ 


56 


0.227 



pressure in a phylogenetically representative set of 23 
species. The gene trees of the obtained partial cichlid 
dlx sequences resulted in various polytomies (Additional 
file 5), probably due to the limited size of some of the 
datasets (minimum of 0.7 kb). Although for each gene 
tree specific branches were observed with relative long 
branches, there is not a particular species or clade that 
has evolved under faster rates of evolution in all of the 
dlx loci examined. Interestingly, three branches have 
relative long branch lengths in multiple topologies: the 
branch towards the Lamprologini (dlx2a, dlx4a and dlxSa), 
C. leptosoma (dlx3b, dlx4a and dbcSd) and C.furcifer {dlxla 
and dlx6a). The relative tree lengths (Additional file 5 
and Table 2) of these gene trees reveal similar results as 
the teleost protein trees, with dlx4b and dlx3a evolving 
fastest and dlxla and dlx6a evolving more slowly. 

Observed signatures of selection pressure in cichlid dlx loci 

To investigate signatures of selection pressure in the dlx 
loci, we performed detailed analyses of the d^lds ratios. 
Maximum likelihood parameter estimations for co, p 0,1,2 
and p and q under different evolutionary models can be 
found in Table 3 for all eight dlx loci. Estimations of o 
under the M0 model suggest that the dlx genes evolved 
under purifying selection with co ranging from 0.0001 
(dlxSa) to 0.457 (dlx2a). A small proportion of sites, 



Table 2 Overall and relative tree lengths of cichlid dlx 



gene trees 


dlx paralog 


Length 
(# base pairs) 


Tree length 
(absolute) 


Tree length 
(relative) 10" 4 


dlxla 


737 


0.036 


0.483 


dlx2a 


1371 


0.094 


0.684 


dlx3a 


666 


0.061 


0.910 


dlx3b 


1972 


0.120 


0.609 


dlx4a 


1166 


0.104 


0.892 


dlx4b 


722 


0.068 


0.937 


dlx5a 


1538 


0.093 


0.607 


dlx6a 


1710 


0.093 


0.542 



0.00001-24.2%, was estimated to have evolved neutrally 
(w = 1) under the Mia model. By using models that allow 
(0 to vary among sites, 0.7-12.3% of sites was detected with 
cj > 1 in dlxla, dlx2a, dlx3a, dlx4b and dlx6a. Overall, 
most sites are estimated to have evolved under purifying 
selection, with highest proportions found in dlx3b, dlx4a 
and dlxSa. 

Likelihood ratio tests of the subsequent model com- 
parisons (Table 4) resulted in the rejection of the null 
models in only the following comparisons per loci: dlxla 
(M8a-M8), dlxla (all four comparisons), dlx3a (M0-M3; 
M8a-M8) and dlx4b (all four comparisons). Positively 
selected sites were detected with the BEB in dlx2a (5 sites), 
dlx3a (1 site) and dlx4b (3 sites; see Table 4, Figure 2). The 
less constraining analyses with the NEB resulted in two 
more putative positively selected sites in dlxla (1) and dlx2a 
(1; Figure 2). Fewer positively selected sites were identified 
by the SLR analyses for dlx2a (position: 36; significance: 
99%), dlx3a (37, 157; 99%, 95%) and dlx4b (145; 99%). 

None of the performed LTRs of the branch-site analyses 
were significant (1 > p > 0.20) indicating that although the 
cj ratios do vary among sites (see above), the <a ratios do 
not vary significantly among lineages. 

Amino acid substitutions and their predicted effect 
on function 

Next, the individual amino acid substitutions were ex- 
amined in more detail. The total protein length and the 
number of amino acid substitutions per locus are shown 
in Table 5 (see also Figure 2 and Additional file 6). A 
relative large number of substitutions was observed in 
dlx2a (13), dlx3a (16) and dlx4b (10), while in dlx5a no 
substitution was found. Most of the amino acid substitu- 
tions are species-specific (i.e., observed in a single species), 
although lineage-specific substitutions were observed for 
the lamprologines (dlx2a, dlx3a, dlx4b), ectodines (dlx2a) 
and haplochromines (dlx2a, dlx3a). None of the observed 
amino acid substitutions have a predicted effect on the 
protein functions (see Table 5), although two substitu- 
tions were observed in the homeobox domain of dlx2a 
(Figure 2). 

Selection regimes on the dlx clusters 

It is known that the paired members of each of the four 
dlx clusters (Additional file 4) are transcribed concur- 
rently [35,42-45]. To characterize if the members of the 
same dlx cluster evolved at similar rates and under simi- 
lar selection regimes, we had a closer inspection of these 
paired genes. First, the teleost dlx protein and cichlid 
gene trees show that overall and relative tree lengths (or 
the rate of evolution) differ between the two genes within 
a cluster. Loci with the highest (dlx3a: 0.583/0.910 and 
dlx4b: 0.864/0.937) or the smallest {dlxla: 0.366/0.483 
and dlx6a: 0.227/0.542) tree lengths are never observed 



Table 3 Site model parameter estimates generated by the CodeML analyses for the eight dlx paralogs 



Locus Parameter estimates under different models 





MO (one ratio) 


Mia (neutral) 


M2a (selection) 


IV13 (discrete) 


M7 (B) 


M8 (fi & b>) 


M8a (B & co s = 1) 


dlx 1a 


to = 0.111 


p 0 = 0.940, to o = 0 


p 0 = 0.980, to„ = 0 


p 0 = 0.640, to o = 0 


p = 0.005, q = 0.049 


p = 0.005, q = 22.651 


p= 1.333, q = 24.671 






p,= 0.060, to, = 1 


p, = 0, oj, = 1 


p, = 0.340, tu, = 0 




p 0 = 0.980 


p 0 = 0.936 








p 2 = 0.020, w 2 = 12.580 


P2 = 0.020, 0)2=12.580 




p, = 0.020,0) = 12.580 


p, = 0.065, oj = 1 


dlx2a 


oj = 0.457 


p 0 = 0.833, oj 0 = 0 


p 0 = 0.879, to o = 0 


p 0 = 0, 0) 0 = 0 


p = 0.005, q = 0.011 


p = 0.005, q = 0.046 


p = 0.005, q= 12.618 






p, = 0.167, w, = 1 


p, = 0.1 10, to, = 1 


p, = 0.986, to, = 0.081 




p 0 = 0.989 


p 0 = 0.833 








p 2 = 0.010, O) 2 = 18.903 


P2 = 0.01 4, 0)2=16.629 




p, = 0.01 1,0) = 18.053 


p, = 0.167, oj= 1 


dlx3a 


oj = 0.320 


p 0 = 0.758, oj 0 = 0 


p„ = 0.921, to o = 0 


p 0 = 0.877, U) 0 = 0 


p = 0.005, q = 0.01 7 


p = 0.01 5, q = 0.088 


p = 0.005, q = 1 3.826 






p, = 0.242, to, = 1 


p, = 0.034, to, = 1 


p, = 0.120, 0), = 2.099 




p 0 = 0.032 


p 0 = 0.758 








p 2 = 0.045, 0)2 = 4.451 


P2 = 0.003, 0)2 = 1 1 .657 




p, = 0.032, 0) = 5.096 


p, = 0.242, oj = 1 


dlx3b 


a) = 0.047 


p 0 = 1, qj 0 = 0.047 


p 0 =1, to 0 = 0.047 


p 0 = 0.243, cu 0 = 0.047 


p = 4.896, q = 99.00 


p = 4.895, q = 99.00 


p = 4.895, q = 99.00 






p, = 1 E-5, to, = 1 


p, = 0, to, = 1 


p, = 0.262, to, = 0.047 




Po= 1 


Po= 1 








P2 = 0, o) 2 = 1 


p 2 = 0.495, tu 2 = 0.047 




p, = 1E-5,to = 1 


p, = 1E-5, to = 1 


dlx4a 


oj = 0.050 


p 0 = 1, oj 0 = 0.050 


p 0 =1, oj 0 = 0.050 


p 0 = 0.216, too = 0.050 


p = 5.267, q = 99.00 


p = 5.266, q = 99.00 


p = 0.1 37, q = 2.305 






P; = 1 E-5, to, = 1 


p, = 0, to, = 1 


p, = 0.384, to, = 0.050 




Po= 1 


Po = 1 








P2 = 0, toj = 1 


P2 = 0.400, u>2 = 0.050 




p, = 1 E-5, oj = 1 


p, = 1 E-5, oj = 1 


dlx4b 


(D = 0.259 


p 0 = 0.872, to o = 0 


p 0 = 0.993, oj 0 = 0.090 


p 0 = 0, oj 0 = 0 


p = 0.005, q = 0.042 


p = 9.871, q = 99.00 


p = 0.005, q = 30.81 7 






p, = 0.128, to, = 1 


p, = 0, oj, = 1 


p, = 0.993, to, = 0.090 




p 0 = 0.993 


p 0 = 0.872 








P2 = 0.007, 0)2 = 12.858 


P2 = 0.007, 0)2=12.858 




p, = 0.007,0) = 12.859 


p, = 0.128, to = 1 


dlxSa 


oj = 1 E-4 


p 0 = 1, to o = 0 


p 0 = 1, aj o = 0 


p 0 = 0.424, aj 0 = 0 


p = 0.005, q = 1 .294 


p = 0.005, q = 99.00 


p = 0.005, q = 44.274 






p, = 1 E-5, to, = 1 


p, = 0, to, = 1 


p, = 0.420, tu, = 0 




Po= 1 


Po= 1 








P2 = 0, (L>2 = 1 


p 2 = 0.156, oj 2 = 0 




p, = lE-5,to = 1 


p, = 1E-5, w = 1 


dlx6a 


a) = 0.029 


p 0 = 0.967, to o = 0 


Po = 0.992, oj 0 = 0 


p 0 = 0.899, u> 0 = 0 


p = 0.066, q = 1 .696 


p = 0.005, q = 2.390 


p = 0.005, q = 1 .672 






p, = 0.033, to, = 1 


p, = 0, to, = 1 


p, = 0.093, oj, = 0 




p 0 = 0.992 


Po = 0.967 








P2 = 0.008, 0)2 = 6.336 


P2 = 0.008, 0)2 = 6.336 




p, = 0.008, O) = 6.336 


p, = 0.033, oj= 1 



Note: P0-2 are the proportions of sites assigned to an oj category or to a beta distribution with p and q as parameters, oj ratios greater than one and their corresponding proportions are depicted in bold. 
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Table 4 Likelihood ratio test (LTR) statistics of site model 
comparisons for dlxla, dlx2a, dlx3a and dlx4b 



Locus 


Test 


LRT (2AZ) 


P 


Selected sites (BEB) 


dlxla 


MO vs M3 


8.416 


0.077 


- 




Mia vs M2a 


3.396 


ns 


- 




M7 vs M8 


3.680 


ns 


- 




M8a vs M8 


5.084 


0.012 


- 


dlxla 


MO vs M3 


70.438 


<0.001 


- 




Mia vs M2a 


39.198 


<0.001 


36A, 48 T, 255/ 




M7 vs M8 


43.093 


<0.001 


36A, 48 T, 195 A, 254 A, 2551 




M8a vs M8 


39.168 


<0.001 


See M7 vs M8 comparison 


dlx3a 


MO vs M3 


1 2.605 


0.013 






Mia vs M2a 


3.858 


ns 






M7 vs M8 


4.258 


ns 






M8a vs M8 


3.872 


0.025 


375 


dlx4b 


MO vs M3 


39.110 


<0.001 






Mia vs M2a 


16.940 


<0.001 


48Q, 135 A, 145 T 




M7 vs M8 


17.367 


<0.001 


48Q, 135 A, 145 T 




M8a vs M8 


16.931 


<0.001 


See M7 vs M8 comparison 



LRT values, p-values and positively selected sites identified by the BEB (p < 0.01 
(in bold) and p < 0.05 (in italic); CodeML) are shown. 



within the same cluster. Furthermore, the mode of selec- 
tion seems to differ between members of the same dlx 
clusters as well. While strong purifying selection was 
observed for dlx3b, dlx4a and dlxSa, their paired cluster 
members dlx4b, dlx3a and dlx6a show sign of elevated 
w-values. A notable exception to this observation is the 
dlxla-dlx2a cluster. For both genes a proportion of sites 
was found with elevated w-values (note that the propor- 
tion is considerably bigger for dlx2a). These observations 
indicate that although clusters are transcribed concur- 
rently, selection seems to act on the individual gene level 
rather than on the level of the dlx gene clusters. Also the 
observed patterns are not in concordance with the two 
groups of homeobox domains that emerged from the 
initial tandem duplication (see [52] and Additional file 4). 

Discussion 

In this work, we present a detailed evolutionary charac- 
terization of the dlx gene repertoire in East African cichlid 
fishes. Previously, Renz et al. [35] studied the embryonic 
expression patterns of dlx genes in cichlids and showed 
that they are expressed in e.g., the developing jaw appar- 
atus and anal fin, tissues that contribute to two putative 
evolutionary innovations: the pharyngeal jaw and the 
egg-spots on the anal fin of the cichlid A. burtoni. Here, 
we study the molecular evolution of dlx genes in a rep- 
resentative set of 23 East African cichlid species. We 
performed comparative phylogenetic analyses and de- 
tailed screens of nonsynonymous-synonymous substitu- 
tion rate ratios to determine the selective pressure acting 



upon these candidate genes for evolutionary novelties in 
cichlid fishes. 

Dlx3a did not get lost in the cichlid lineage 

Our phylogenetic analyses of dfx proteins extends previous 
analyses (e.g., [35]) by the inclusion of cod [72] and four 
different cichlid species (i.e., O. niloticus, N. brichardi, 
A. burtoni and P. nyererei; BROAD Institute). Although 
our results agree with most of the available hypotheses 
on the evolutionary loss of dlx genes in specific teleost 
lineages (i.e., dlx3a in zebrafish and dlx4a in medaka), 
we did detect dlx3a in cichlids and thus refute the 
cichlid-specific gene loss hypothesis of dlx3a put forward 
by Renz et al. [35]. Not only were we able to locate this 
gene in all four cichlid genomes examined (Additional 
file 1), we also gathered partial gene sequences for this 
locus in all 23 cichlid species included (Additional files 4 
and 6). Furthermore, in-house tblastx searches of this 
newly identified paralog against preliminary cichlid EST 
libraries (BROAD Institute, unpublished data) resulted 
in multiple hits, providing proof of its expression in - at 
least - Astatotilapia burtoni, Oreochromis niloticus and 
Metriaclima zebra. 

Selection on dlx paralogs in relation to gene 
duplication events 

Gene-wide estimates of the d N /d s ratios indicate that all 
loci evolved under purifying selection (cj < 1), indicating 
strong selection against deleterious mutations, commonly 
observed in functional proteins. Additional analyses of 
individual codons indicate that the sequenced regions of 
dlx3b, dlx4a and dlxSa evolved under purifying selec- 
tion, while positive selection acting on specific codons 
was detected for a small proportion of sites (i.e., up to 
12%) for dlxla, dlx2a, dlx3a and dlx4b (i.e., a smaller 
number of positively selected sites was found with more 
stringent SLR analyses for dlx2a, dlx3a and dlx4b). 
Plausible reasons for the excess of nonsynonymous muta- 
tions in these loci are either lowered functional constraints 
or directional selection, as Sumiyama and colleagues sug- 
gested for Dlx7 in mouse [73]. Different modes of selec- 
tion are thus found to have acted on the dlx paralogs in 
cichlids after the genome duplication events. 

Differential selection after gen(om)e duplication is a 
commonly observed phenomenon and is associated with 
the fate of the gene duplicates i.e., non-, sub- or neofunc- 
tionalization. Sub- and neofunctionalization are adaptive 
processes by which either spatial or temporal partitioning 
of the ancestral function or the evolution of complete new 
functions take place [5,74-76]. While ancestral functions 
can be maintained by retaining the protein sequences and 
preventing deleterious mutations through purifying selec- 
tion, relaxed selection on the other duplicate can lead to 
the introduction of mutations and subsequent divergence 



Diepeveen et al, BMC Evolutionary Biology 2013, 13:153 
http://www.biomedcentral.com/1471 -21 48/1 3/1 53 



Page 8 of 12 



nf : llllllllllllllllllllllllllimllllllllllll C 



IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIM 



Pred: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA: GEPPEEMTKTTIPESLNSPASGKTVFMEFGPPSQQMSPSS 



Pred : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA: CSISASKDSPTMPESSSTDMGFYSGQSALHGSQDFYlAQQ 



llllllllllllllllllllllllllllllllllllllll 



cccccccccccccccccccccccccccccccccccccccc 

MSHGHYPMHCLHSAGHTQHDSYSPASSFPRSLGYPYVNSV 



Conf ]||||||||||||||||||||||||||||||||||||||||E 

Pred: 

Pred: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA: PYSAQHMNPYAYHHYSLNGMCPGCAYPVGKAEYPYPHAAY 



o f iiimiimiimiimimimmmii^if 

red : 

red: CCCCCCCCCCGCCCCCCCCCCCCGCCCCCCCCCCCCGCCG 
AA : GSHSTSPYLSTVQTYQNSSGLTQTRLE0AAPETEKKTWE 

90 100 110 12C 

onf: 



^■■■■□^□^□□■lllllllllllllllllllllllllE 



90 100 110 120 

conf iiiiinniiiimimiiiiimiimmiimi 



iiiiiimiiiiiiiimmmiiiiiiimiiii 



Pred : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA: Q I TSSN YH SL HKS QE S P T LP V5 T ATD S S Y YNGQQPfcc AG 



]||||||||||||||||||||||||||||||IHIIIIIII 



Pred: 

AA: SPFGQLgIyQYHSSATSSVPYNAKSYDLGFNSSYGTYGSY 
** ' * 1 * * 1 



red: CCCCCCCCO 
AA: PPSPAVWENNNSSNGNNNSSNHSSGSGNNPNVSNGSfcNS 

130 * 140 150 * * *160 

onf llllllllllllllllllllllllllllllllllllllll 



Pred: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA : AMLDPTSRGQVPYQPPVDSPPAYMGEYTHQNWYQQQGAH L, 
* * * ** * ' 

170 180 190 200 



llllllllllllllllliilllE 



red: CCCCCCCCCCCCCCCCCCCCCC 
AA: GLSQSGQVHHAPPTAPSATQSM 



nf jlllllllllllllllll^r, ullllllllllllllllll D 



-c- 



-n 



Conf : IlllMllllllllllllIlm^llllllllllllllll 

Pred: strand Pred: predicted secondary structur e 

Pred : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA; MMSVGFIPDSLNGSDPSKSAFLEFGHGHSTHQQHSSGLSH 



illlllllll: rjllli^llllllillllli^^iir 

Pred: 



Prod: - i [■ ■ '■ ■ :cc:lhhhh:1i:ccccccchhhcccc 

AA : Q L AALQRRF QKTQYLALPE RAE L AA SLGLTQTQVKIWFQN 



m jlll^^^lllllllllllllllllllllllllllllll 



red : ■ CCCCCCCC 

AA: VYPVHGL^AGHSQHESPFPGTASYGRSLGYAYPGTVNTH 



imiimmiiiimmiimiimiimiiiiE 



-z 



Pred: IHKHHKHHiiHCCCCCCCCCCC' :G 

AA: ^^^^^pKSGEIPPEQHVASSESPPCTSPPTT^IDFPQ 

170 180 190 20C 

Conf : jllllllllllllllllllllllllllllllllllllllll 

Pred: 

Pred: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
AA : TQRMNNVSSSLPQSBSPPNTTAPSSFLANYSWYSTTNSTT 

210 220 230 24C 

Conf: ]|||||||||||||||E 



o 1 ODllllllllllllllllllll 

red: 

red: CCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCC 
AA: PPSAYMPYQX^^^^MkHGSGSEGEHLHGTSSISI-CS? 
* 

90 100 110 120 

onf: lllllMllllllllllllllllllllllllllllllllll 

red : 

red: CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
A A : I I H II :YXI\H' . I Y : ' M : : . 1 1 -■ I l . 1 0 !A 

*'***' 

130 140 150 160 
Conf: 1IIIIIIE 



:..,}■ CCCCCCCC' 
AA: HLQPPLVQHHHNSB 
** 

250 



I = homeobox domain. | 
x C„,:l: L 3 ]||[ 



i positive f 

^nfidence a 



Elected site 



red: CCCCCCC 
AA: MHRPQMM 



Figure 2 Secondary structure and positively selected sites for four partially sequenced Astatotilapia burtoni Dlx proteins. Secondary structure 
predictions were obtained from the PSIPRED server (http//bioinf.cs.ucl.ac.uk / psipred/). Positively selected sites identified by the site model analyses 
(CodeML) and the SLR analyses are highlighted in red (BEB and/or SLR) or orange (NEB) boxes. (A) Dlxla. (B) Dlx2a. (C) Dlx3a. (D) Dlx4b. 



[5,75,76]. Most of these changes are deleterious and are 
followed by the loss of the gene over time (i.e., nonfunc- 
tionalization). On rare occasions the mutations can lead 
to an altered function of the protein (i.e., neofunctional- 
ization; change within the protein) or altered expression 
pattern (subfunctionalization; change in regulatory regions), 
which can be characterized by elevated cj values and the 
maintenance of the mutations results in divergence of the 
two duplicates. 



Many studies have focused on duplicated genes in rela- 
tion to divergence of duplicates (see e.g., [77-80] and refer- 
ences therein). An interesting case of subfunctionalization 
was described in leaf-eating Colobine monkeys, in which 
the pancreatic ribonuclease gene (RNASE1), necessary 
to digest its specialized diet, was duplicated [81,82]. Al- 
though the two gene-products are used in the same 
process (i.e., digestion of bacterial RNA), the duplicate 
gene shows many substitutions, while the ancestral locus 
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Table 5 Amino acid substitution and their predicted 
effect on function for the eight cichlid dlx loci 



Locus 


L 


S 


E 


dlx la 


144 


1 


Tolerated 


dlx2a 


255 


13 


Tolerated 


dlx3a 


222 


16 


Tolerated 


dlx3b 


160 


4 


Tolerated 


dlx4a 


96 


1 


Tolerated 


dlx4b 


167 


10 


Tolerated 


dlx5a 


271 


0 




dlx6a 


236 


2 


Tolerated 



Length of obtained protein sequence (L), the number of amino acid 
substitutions (S) and the predicted effect of the substitutions on the function 
generated by the SIFT analyses (E) are shown. 



did not change [81]. Similar patterns of heterogeneity 
in amino acid substitutions or differential selection were 
also observed by Dermitzakis and Clark [83] between 
duplicates of several developmental gene families (e.g., 
Notch, Bmp and Hox9) in mouse and human. Interest- 
ingly, differential selection regimes acting on paralogs 
were also found in the murine Dlx3-Dlxl cluster, with 
Dlx7 evolving more rapidly than Dlx3 [73]. Our results 
of differential selection acting on the cichlid dlx paralogs 
are thus comparable to previously studied cases of du- 
plicated genes. We even detect a similar pattern as 
Sumiyama et al. [73], with dlx4b evolving more rapidly 
than dlx3b (i.e., relative tree length 0.937 vs 0.609). 

The adaptive protein evolution as observed in dlxla, 
dlx2a, dlx3a and dlx4b together with the evolutionary 
history of the gene family, could thus be a sign of pos- 
sible new or altered functions of these dlx paralogs in 
cichlids. Although we did not observe amino acid substi- 
tutions with predicted apparent effect on the protein 
function in our partial sequences, other mechanisms, such 
as as-regulatory evolution might have altered the expres- 
sion patterns after gene duplication. Gene expression 
analyses in cichlid and zebrafish indicate that clusters 
are often transcribed concurrently and that the dlx du- 
plicates exhibit overlapping expression patterns in par- 
ticular during the development of brain and pharyngeal 
arches [35,44,45]. This co-expression of the dlx clusters 
is controlled through intergenic c/s-regulatory regions 
[35,42,43]. While mutations in these regions are expected 
to affect the expression of both paralogs, changes in the 
coding regions of the dlx loci are likely to affect the 
individual dlx locus' function, which could lead to 
neofunctionalization. 

Selection pressure on dlx paralogs in relation to 
evolutionary innovations 

We found an interesting pattern comparing our d^/ds 
results with the expression patterns found by Renz et al. 



[35] in relation to evolutionary novelties in cichlids. In 
the developing pharyngeal teeth and the anal fin dlx3b, 
dlx4a (not in anal fin) and dlxSa, the exact loci for 
which we found strong patterns of purifying selection, 
are expressed at high levels. Although this observation 
seems to contradict other cases in which candidate 
genes showed accelerated rates of protein evolution (see 
[31,53,84]), they do not stand alone (see e.g., [36]). It 
has been shown that minor changes in the complex gen- 
etic pathways underlying the development of morpho- 
logical structures can lead to the evolution of novelties 
(see e.g., [85]). Furthermore, many cases of morpho- 
logical adaptation are driven by cw-regulatory evolution 
(reviewed in [86]). Several intergenic ci's-regulatory ele- 
ments have been identified in the dlx clusters in A. burtoni 
by Renz et al. [35], but the functional characterization 
in cichlids is yet to be performed. It is thus possible that 
only a small fraction of genes involved in the evolution- 
ary novelties in cichlids show signs of adaptive evolution 
and that the three dlx loci were co-opted for their an- 
cestral functions. 

According to Renz et al. [35], the five dlx genes for 
which we found signatures of positive selection, are either 
not expressed at all or at low levels during pharyngeal 
teeth and anal fin development in the cichlid A. burtoni. 
Low levels of gene expression were observed for dlx2a in 
the developing pharyngeal teeth in cichlids [35], while 
higher dlx2a expression levels were observed in other 
teleost species [33,44,45]. Dlx4b and dlx6a expression 
has previously been shown in the developing pharyngeal 
teeth of zebrafish and/or medaka [44,45], but has not been 
observed in cichlids (yet). Furthermore, multiple dlx genes, 
including loci with signatures of positive selection, appear 
to be expressed in the developing anal fin tissue at time 
points coinciding with egg-spot development in A. burtoni 
(E. Santos, personal communications). Therefore, it is likely 
that several dlx paralogs, for which we found signs of 
positive selection, are involved in the development of 
evolutionary innovations in cichlids, in contrast to the 
initial findings of Renz et al. [35]. Future detailed and 
extended functional analyses should be conducted to 
elucidate their role in the development of these evolu- 
tionary important traits in cichlid fishes. 

Conclusions 

In this study, we provide an in depth molecular evolu- 
tionary analysis of the dlx gene repertoire in teleost 
fishes. We located and generated partial sequences for 
dlx3a in 23 East African cichlid species, refuting the hy- 
pothesis of Renz et al. [35] that dlx3a got lost in the 
cichlid lineage. Phylogenetic analyses of the teleost dlx 
gene repertoire show that substantial differences exist 
in the rate of evolution among teleost dlx paralogs. In 
addition, analyses of the nonsynonymous-synonymous 
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substitution rates of the cichlid dlx paralogs revealed 
strong differences in the selection pressure acting upon 
dlx paralogs and cluster members. Although differen- 
tial selection pressure after gene duplication is a puta- 
tive sign of new or altered functions, we observed a link 
between the dlx loci under strong purifying selection, in 
particular, and high expression levels in two cichlids' nov- 
elties; the pharyngeal jaw and anal fin. This indicates that 
other mechanisms than adaptive protein evolution are 
likely to be involved in the co-option of these genes. Fur- 
thermore, several (preliminary) studies found that at least 
three other dlx paralogs, for which we found signs of 
positive selection, are actually expressed in the developing 
pharyngeal teeth and/or haplochromine anal fin. Hence, 
the dlx paralogs appear as candidate genes for the devel- 
opment of evolutionary innovations in cichlids, although 
further functional analyses should elucidate the role of 
positive selection therein. 

Availability of supporting data 

The datasets supporting the results of this article are pub- 
licly available in the GenBank repository under accession 
numbers: KC285366-KC285546 and in the treeBASE re- 
pository under submission number 14433, http://purl.org/ 
phylo/treebase/phylows/study/TB2:S14433. 

Additional files 



resolved tree, (f) Dlx4b (722 bp; TPM3uf). Mostly polytomous relationships 
between species, except the Limnochromini and most members of the 
Lamprologines. (g) Dlx5a (1 538 bp; TIM2 + G). Basal polytomy divides 
ingroup species except 6. permaxillaris, in two big clades. (h) Dlx6a 
(1710 bp;TIM3 + G). Limnochromines, Lamprologines and 
Haplochromines recovered as monophyletic clades, although the 
relationships between lineages are largely polytomous. 

Additional file 6: Four partially sequenced cichlid Dlx proteins. 

Depicted are the amino acid sequences of Astatotilapia burtoni (a, c, d) 
and Ctenochromis horei (b). Secondary structure predictions were obtained 
from the PSIPRED server (http://bioinf.cs.ucl.ac.uk/psipred/). (a) Dlx3b. 
(b) Dlx4a. (c) Dlx5a. (d) Dlx6a. 
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Additional file 1: Accession numbers and/or genomic location of 
the teleost dlx sequences. 

Additional file 2: Specimen information and GenBank 
Accession numbers. 

Additional file 3: Primer information and primer sequences. 

Additional file 4: Protein comparison of the teleost dlx homeobox 
domains. Depicted are the amino acid sequences of the homeobox 
domains for each of the four teleost clusters: dlxla-dlx2a, dlx4a-dlx3a, 
dlx4b-dlx3b and dlx6a-dlx5a in comparison with the single Dll homeobox 
sequence (here depicted in duplo) of Drosophila melanogaster. 
Sequences can be divided in two groups; dlxla, dlx4a, dlx4b and 
dlx6a versus dlx2a, dlx3a, dlx3b and dlx5a. The two sixty amino acid 
long homeobox domains of each cluster are depicted in separate 
boxes. The top graph displays the mean pairwise identity of all 
sequences (i.e., green = 100% identity and brown > 30% identity). 
Numbers represent the amino acid position within the homeobox. 

Additional file 5: Maximum likelihood gene trees based on 23 
cichlid species for the eight dlx loci. Bootstrap values (PAUP*) and 
Bayesian posterior probabilities (MrBayes) above 50% are shown 
respectively above and below the branches. A color key for the ten 
studied cichlid lineages is given in the box below the figure, (a) Dlxla 
(737 base pairs (bp); TPM3uf model). Two major polytomies were 
recovered. The lamprologines cluster together with the 
Boulengerochromini, Bathybatini and the Cyphotilapiini. A. burtoni is 
found at the base with O. tanganicae. (b) Dix2a (1371 bp; HKY + I model). 
Polytomous tree with all members of the lineages Lamprologines, 
Ectodines, Haplochromines and Limnochromines recovered as 
monophyletic clades. (c) Dlx3a (666 bp; HKY model). Polytomous tree, 
with only the Lamprologines recovered as monophyletic clade. (d) Dix4a 
(1 166 bp; TPM3uf + I + G). Polytomous relationships were observed 
between multiple lineages, although most lineages are monophyletic 
except the Haplochromines (e) Dlx3b (1972 bp; GTR + I + G). Moderately 
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